library(ape)

# Func to get height to root
root.height <- function (x) {
	n <- length(x$tip.label)
	n.node <- x$Nnode
	N <- n + n.node
	rootnode <- (1:N)[!(1:N) %in% x$edge[,2]]
	d <- dist.nodes(x)[rootnode, 1:n]
	res <- max(d)
	res
}

true.trees <- read.nexus("exampleLOH_true_trees.trees")
n <- length(true.trees)

# Calc max root height
# Normalize branch lengths to be in units of years
max.rh <- 0
for (i in 1:n){
	true.trees[i]$tree$edge.length <- true.trees[i]$tree$edge.length/300.0
	rh <- root.height(true.trees[i]$tree)
	max.rh <- max(rh, max.rh)
}
max.rh <- max.rh+1

par(mfrow=c(1,1), mai=c(0.2, 0.2, 0.2, 0.2), mar=c(2,2,2,2), new=FALSE)
for (i in 1:n){
 color <- sample(rainbow(20),1)
 plot(x.lim=c(0,max.rh),y.lim=c(0,30),ladderize(beast.trees[i]$tree),edge.color=color,direction="l",type="c",show.tip.label=FALSE,node.pos=1)
 axisPhylo()
 par(new=TRUE)
}
plot(x.lim=c(0,max.rh),y.lim=c(0,30),ladderize(true.tree),edge.color="black",edge.width=3,direction="l",type="c",show.tip.label=FALSE)
axisPhylo()
box(which="plot",lty="solid")



